function [R,E]=mcf_id(A)
% Cholesky factorization with added multiple of the identity;
% based on Algorithm 3.3 in Numerical optimization (NW 2) p51
% INPUT:
%       A - the matrix to factorize;
% OUTPUT:
%       R - the triangular factor matrix;
%       E - the error matrix;
% REVISION:
%       jqin, 01/feb/2011, created;
n=numel(A(:,1));
R=spalloc(n,n,(1+n)*n/2);
E=spalloc(n,n,n);
for k=1:n
    if k~=n
        R(k,k)=sqrt(abs(A(k,k)));
        E(k,k)=R(k,k)^2-A(k,k);
        for j=k+1:n
            R(k,j)=A(k,j)/R(k,k);
            for i=k+1:j
                A(i,j)=A(i,j)-R(k,j)*R(k,i);
            end
        end
    else
        R(k,k)=sqrt(abs(A(k,k)));
        E(k,k)=R(k,k)^2-A(k,k);
    end
end
